function    simulate( sim, endTime )
% Get system status before cutting time via simulation
% 
% If user chooses to reuse older file and the file exists, exit directly.
global dataDir dataID dataSim postSimFunc guessSimFunc
dataSim = strcat(dataDir, dataID, '.mat');

global Y3 Y3fun L_Y3 U_Y3 p_Y3 uEXd uGall
postSimFunc = @postSimulate;
guessSimFunc = @guessSimulate;

global NTswitch sbspara                         % by case$x.m
t = 0; h = sbspara(1);
nN1 = fix(NTswitch(5)/h) +1; nN2 = fix(NTswitch(6)/h) +1; nN3 = fix(endTime/h) +1;
NTswitch(5) = (nN1 -1) *h; NTswitch(6) = (nN2 -1) *h;

if sim == 0 && exist(dataSim, 'file') == 2
    return
end

%% Simulate using Runge-Kutta (RK-4)
% IVP: Dy = f(t, y); y(0) = y0;
%   y(n+1) = y(n) + h/6 * (k1 + 2*k2 + 2*k3 + k4); t(n+1) = t(n) + h
%   where: k1 = f(t(n), y(n));        [at beginning]
%          k2 = f(t(n)+h/2, y(n)+h/2*k1); k3 = f(t(n)+h/2, y(n)+h/2*k2); [at midpoint]
%          k4 = f(t(n)+h, y_k4(n)); [at endpoint ][y_k4(n) = y(n) + h*k3]
global nB nG idx_gx idx_gy iG_cut               % by initial.m
global delta omega omgs Pm EXd TJ TJsum         % by initial.m
global Y1 Y2 I V                                % by Ymatrix.m
global dlt_cut omg_cut V_cut

% matrix decomposition
[L_Y1, U_Y1, p_Y1] = lu(Y1, 'vector');          % L*U = P*A; Ax = b --> x = U\(L\b(p))
[L_Y2, U_Y2, p_Y2] = lu(Y2, 'vector');          % "P --> 'vector' p" is faster
% parameters
tRec = zeros(nN3, 1); deltaRec = zeros(nN3, nG); coiRec = zeros(nN3, 1);
I = zeros(2*nB, 1); V = zeros(2*nB, 1); uEXd = EXd;

%% Phase1: 0 --> fault clear; Phase2: fault clear --> start shedding
simPhase = 1;
for iT = 1 : nN3                                % total discreted elements
    tRec(iT) = t;   t = t + h;
    deltaRec(iT, :) = delta *180 /pi;           % Record history
    coiRec(iT) = sum(delta .* TJ) / TJsum *180 /pi;
    
    if iT == nN1, simPhase = 2;                 % at fault clear time
    elseif iT == nN2, dlt_cut = delta; omg_cut = omega; V_cut = V;
    elseif iT == nN3, break; 
    end
    RungeKutta(simPhase);
end

%% Store result
save(dataSim, 'tRec', 'deltaRec', 'coiRec', ...
    'dlt_cut', 'omg_cut', 'V_cut');

%% functions
    % ODE (delta, omega)      [EXd = E'/Xd']
    function [Ddlt, Domg] = diffEqua(dlt, omg, phase)
        I(idx_gx) = uEXd.*sin(dlt); I(idx_gy) = -uEXd.*cos(dlt); 
        if phase == 1
            V = U_Y1\(L_Y1\I(p_Y1));            % V = Y1 \ I;
        elseif phase == 2
            V = U_Y2\(L_Y2\I(p_Y2));            % V = Y2 \ I;
        elseif phase == 3
            V = U_Y3\(L_Y3\I(p_Y3));            % V = Y3 \ I;
        end
        Pe = EXd.*(V(idx_gx).*sin(dlt)-V(idx_gy).*cos(dlt));
        Ddlt = omgs*(omg - 1);
        Domg = (Pm - Pe)./TJ;
    end

    function RungeKutta (phase)
        % k1 k2 k3 k4
        [k1_dlt, k1_omg] = diffEqua(delta, omega, phase);
        dlt_k2 = delta + h/2 * k1_dlt; omg_k2 = omega + h/2 * k1_omg;
        [k2_dlt, k2_omg] = diffEqua(dlt_k2, omg_k2, phase);
        dlt_k3 = delta + h/2 * k2_dlt; omg_k3 = omega + h/2 * k2_omg;
        [k3_dlt, k3_omg] = diffEqua(dlt_k3, omg_k3, phase);
        dlt_k4 = delta + h * k3_dlt; omg_k4 = omega + h * k3_omg;
        [k4_dlt, k4_omg] = diffEqua(dlt_k4, omg_k4, phase);
        % finish RK-4
        delta = delta + h/6 * (k1_dlt + 2* k2_dlt + 2* k3_dlt + k4_dlt);
        omega = omega + h/6 * (k1_omg + 2* k2_omg + 2* k3_omg + k4_omg);
    end

    function [tPRec, dltPRec, coiPRec, omgPRec] = postSimulate (uG, uL)
        nPRec = nN3 - nN2; t = nN2 * h;
        tPRec = zeros(nPRec, 1); 
        dltPRec = zeros(nPRec, nG); omgPRec = dltPRec;
        coiPRec = zeros(nPRec, 1);
        feval(Y3fun, uG, uL);                           % update Y3 via solution
        [L_Y3, U_Y3, p_Y3] = lu(Y3, 'vector');          % "P --> 'vector' p" is faster
        delta = dlt_cut; omega = omg_cut;               %V = V_cut;
        uGall(iG_cut) = uG; uEXd = (1 - uGall) .* EXd;  % parameters after cutting
        %uTJ = (1 - uGall) .* TJ; uTJsum = sum(uTJ);
        for iPT = 1 : nPRec
            RungeKutta(3);
            tPRec(iPT) = t;   t = t + h;
            dltPRec(iPT, :) = delta *180 /pi;            % Record history
            omgPRec(iPT, :) = omega;
            coiPRec(iPT) = sum(delta .* TJ) / TJsum *180 /pi;
        end
    end

    function [dltGRec, omgGRec, VGRec, coiGRec] = ...
            guessSimulate (tG, uG0, uL0)
        nStep = nN3 - nN2; t = nN2 * h;
        nC1 = length(tG); nC = nC1 - 1;
        dltGRec = zeros(nC1, nG); omgGRec = dltGRec;
        VGRec = zeros(nC, 2*nB); coiGRec = zeros(nC, 1);
        feval(Y3fun, uG0, uL0);                           % update Y3 via solution
        [L_Y3, U_Y3, p_Y3] = lu(Y3, 'vector');          % "P --> 'vector' p" is faster
        delta = dlt_cut; omega = omg_cut;               %V = V_cut;
        uGall(iG_cut) = uG0; uEXd = (1 - uGall) .* EXd;  % parameters after cutting
        %uTJ = (1 - uGall) .* TJ; uTJsum = sum(uTJ);
        iC = 1; dltGRec(iC,:)=delta; omgGRec(iC,:)=omega;
        for iPT = 1 : nStep
            RungeKutta(3); t = t + h;
            if iC+1 > length(tG), break; end
            if t + 1e-5 > tG(iC+1)
               
                    dltGRec(iC+1,:) = delta;
                    omgGRec(iC+1,:) = omega;
                    VGRec(iC,:) = V; 
                    coiGRec(iC) = sum(delta .* TJ) / TJsum;
          
                iC = iC +1;
            end
        end
    end
end